Hamiltonian adaptive resolution simulation for molecular liquids 
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Abstract 

Adaptive resolution schemes allow the simulation of a molecular fluid treating simultaneously 
different subregions of the system at different levels of resolution. In this work we present a new 
scheme formulated in terms of a global Hamiltonian. Within this approach equilibrium states 
corresponding to well defined statistical ensembles can be generated making use of all standard 
Molecular Dynamics or Monte Carlo methods. Models at different resolutions can thus be coupled, 
and thermodynamic equilibrium can be modulated keeping each region at desired pressure or 
density without disrupting the Hamiltonian framework. 
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Complex molecular fluids and soft matter typically display inherently multiscale phe- 
nomena and properties. To handle this problem general strategies have been developed, 
which can be classified as either sequential or simultaneous. In the former class of meth- 
ods, coarse-grained models are (usually) developed from microscopic input l|-l3|; systems 
are then simulated separately at different levels of resolution. In the latter class, which we 
pursue here, systems are treated within a single simulation on different levels of resolution. 
A small, well defined region of space is kept at a higher level of detail, while the surrounding 
can be treated on a coarser, computationally more efficient level. 

This idea has been successfully employed, for example, to investigate crack propagation 
in hard matter ^-8] and in mixed quantum mechanics/molecular mechanics (QM/MM) sim- 
ulations, where particles are assigned statically to either the QM or the MM region 9Nl3|. 
For soft matter and liquids, inherent fluctuations and particle diffusion require a setup where 
molecules, or parts of them, can cross boundaries between areas at different resolution, while 
maintaining the overall thermodynamic equilibrium. Scale-bridging methods have been de- 
veloped in various fashions to couple all atom (AA) and coarse-grained (CG) models fl^ . 



and even particle-based models to the continuum 



ISMlTj]. To our knowledge, to date the 



only energy-conserving mixed-resolution approach is the 'Adaptive Partitioning of the La- 



grangian' method by Heyden and Truhlar 



19| . In this method, a combinatoric sum of all 



possible AA and CG interactions between molecules in different resolution regions is taken 
into account. The practical viability of this approach is limited by its intrinsic combinatoric 
complexity, and by the fact that the resulting equations of motion are not amenable to a 
standard symplectic integrator (e.g. velocity Verlet of leap frog), so that an ad hoc, more 
complicated one had to be developed. 

With this idea of mixed resolution in mind the AdResS (Adaptive Resolution Scheme) 
method was developed, in which one can dynamically couple specific regions of a simulation 



box at different leve 
rium between them 



resolution, while maintaining the correct thermodynamic equilib- 
26|. The particles move from one region to the other through a 



s oi reson 
3, 3-26 



hybrid resolution zone (Fig. [T]): in this region the resolution switch is defined by a transition 
function A(x), smoothly changing the interactions from an atomistic description, A = 1, to 
a coarser one, A = 0, which typically contains a considerably smaller number of degrees 
of freedom (DOF's) per molecule. AdResS is based on the requirement that molecules 
interact through pairwise forces, and Newton's third law is strictly satisfied in the whole 
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FIG. 1. Illustration of the simulation setup for adaptive resolution simulations. Molecules freely 
move from an atomistic region, AA, A = 1, through a transition zone H to the coarse-grained 
region, CG, A = 0. A(R), R being the position of the center of mass of the molecules, is a smooth 
transition function that interpolates between the AA and the CG region. 



simulation box by construction. These requirements lead to a force interpolation scheme 
between molecules, Fq,/^ = A(Rq)A(R/3)F;^^ + [1 - K^a)H^i3)]Fa^ , where the force 
between centers of mass of molecule a and /3 consists of an atomistic, A(Ra)A(Rfl)F^*^ and 
a coarse-grained part, [1 — A(Rq)A(R/3)]F^^. Yet, it was formally demonstrated |27| that a 
Hamiltonian compatible with this force interpolation scheme can not exist. 

The method is nonetheless robust, since it allows us to define temperature, pressure and 
density everywhere, and the introduction of a thermodynamic force 



in the transition 



zone paved the way to open system MD simulations [20|, |29|]. Despite the success of the 
force-interpolation based AdResS method, though, the lack of a Hamiltonian description in 
the transition region is a drawback: it hampers a general statistical theory for the whole 
setup, limits the choice of the simulation ensemble and prevents Monte Carlo simulations. 
Moreover, in the transition region the system has to be stabilized by a local thermostat that 
removes excess heat thus keeping the system in a state of dynamical equilibrium 
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31|. In this Letter we propose a new resolution- interpolation and coupling 



concept, if- AdResS, which is formulated in terms of a general Hamiltonian H for the whole 
system. Furthermore, we develop an analogy to the Kirkwood 32] coupling parameter 
scheme where we relate the variation of the thermodynamic properties through the transition 
zone in H-AdResS to the integration over A in homogeneous systems. We demonstrate 
our approach on a prototypical mixed AA/CG system, showing that the existence of a 
global Hamiltonian makes it possible to perform microcanonical (NVE) adaptive resolution 



simulations. 



Let us consider a system composed of molecules. [33| (labeled by greek indices), each 
having n atoms (labeled by latin indices) The resulting M = nN atoms interact via general 
intramolecular potentials, and pairwise intermolecular potentials. The Hamiltonian of this 
system can be written as: 



TV n 9 N 



Q=l i=l a=l 
^ N n 



I3^a ij 

ytnt indicates the intramolecular interaction, for which we do not need to make any assump- 
tion, pai, fnai, and Vai are the momentum, mass, and position, respectively, of atom i of 
molecule a. We now consider a CG pair potential V^*^ = V^*"*^(Rq, — R^) that depends 
on the center of mass (CoM) positions R of the molecules a and /3; the total CG potential 
energy of molecule a is given by V^'^ = J2i3^a K^*^/^- In analogy to Kirkwood's thermody- 



namic integration (TI) method to compute free energy differences 32|, we define a 'mixed 
resolution' Hamiltonian H: 

H= ^ (2) 

OLi a 

where coupling parameters = A(Ra) were introduced, which depend on the position 
of the molecules' centers of mass. The local resolution A(R) varies between 1 (completely 
AA system) and (completely CG system); intermediate A values define a hybrid region, 
interfacing the atomistic and the coarse-grained ones (as in Fig. [T]). According to the 
Hamiltonian in Eq. [21 for molecules interacting with a mixed-resolution (having A G (0, 1)) 
both AA and CG total potential energies are calculated and weighted according to their 



own resolution A. The atomistic DOF's are retained everywhere |20[ and their dynamics is 
seamlessly evolved, allowing coarse-grained molecules close to the hybrid region to interact 
also at the atomistic level. Being defined in terms of a Hamiltonian, this hybrid-resolution 
scheme conserves the total energy in a microcanonical simulation, as it is numerically verified 



34( 1 . Furthermore, as different regions exchange particles and energy the resulting stationary 



state is an equilibrium state. 



The force derived from H (Eq. [2]) has the form: 



(3) 




— — 2^ * + ^ — * "^1/5 



where F^^^^^- (resp. F^j^) is the AA (resp. CG) force acting on atom i of molecule a due 
to the interaction with molecule /3. The distribution of CoM forces onto the atoms is de- 

The first term of Eq. [3] contains a weighted sum of pairwise forces 



scribed in Ref. 



and is antisymmetric by exchange of the molecules' labels; this term therefore complies 
with Newton's 3rd law everywhere, and is analogous to the AdResS force interpolation. 
The second term, F™*, is the force exerted on atom ai by the other atoms in the same 
molecule and does not contribute to the force balance between molecules. The third term, 
F^^ = — \y^^ — V^'^~\ VQiA(Ro), introduces in the hybrid region a drift force which violates 
Newton's 3rd law and momentum conservation. F^^ plays the role of an external force induc- 
ing, in general, pressure and density inhomogeneities in the system as it reaches equilibrium 
(no temperature gradients are present because energy is conserved and freely flows between 
the two subdomains). In particular, the drift force is balanced by a hydrostatic pressure 



gradient across the hybrid region given by Vj) = p(F'^^) 35 1. 

We validated our approach on the same model system as in (20| and illustrated in Fig^ 
A detailed description of these simulations is given in the supplementary information |34|. 
The AA system consists of tetrahedral molecules, each composed by four atoms of unit 
mass connected by anharmonic springs. The atomistic interaction between molecules is 
given by a purely repulsive Weeks Chandler Andersen (WCA) potential, while the CG 
potential was obtained via Iterative Boltzmann Inversion (IBI) 36|. In contrast to the 
original AdResS scheme, we could perform these adaptive resolution simulations in the 
microcanonical ensemble, achieving conservation of the total energy (Fig. SI). The resulting 
density profiles are fiat in both the AA and the CG regions, and within 1% of the reference 
values (Fig. S2). In addition, the structure of the fluid in the AA region is unchanged 
compared to a purely AA simulation (Fig. S2). 

Such seamless coupling stems from the good matching between the thermodynamic prop- 
erties of the AA potential and the CG potential. In general, however, the employed CG 



potentials are only approximations to the the exact many-body CG potential 37 



. As 



a consequence, there is usually a thermodynamic mismatch between the AA and the CG 
systems, which can have different chemical potentials and equations of state jsl- If compen- 
sations are introduced as time-independent functions of the position of the molecules, then 
the modified Hamiltonian H is expressed as 



TV 



(4) 



and conserves the energy. The compensation terms change the drift force to 

d AH{X) 



AA 



v: 



CG 



V«A(R«) 



(5) 



A=A{R„)_ 

In the following, we relate suitable compensations to the Kirkwood's TI scheme for the 
free energy difference AF{X) between a hybrid system with a position-independent coupling 
parameter A < 1 and a coarse-grained system (A = 0) at the reference concentration p*: 

AF(A) ' 



N 



1 r ,wiAF(AO 
d H{X') 



A 



1 

1 

NJo 



dX' 



dX' 
dX' { [y^^ - 



(6) 



Consider first a situation, where we wish to embed the AA region in a CG region with 
identical molecular (Virial plus kinetic) pressure. |39| To avoid the buildup of a hydrostatic 



pressure gradient 



35 | across the hybrid region, we need to assure that Vp = p{F'^^) = or 
d AH{X) 



dX 



(7) 

A=A(R<,) 

If we replace the local average at each given A = A(Rq) by the corresponding value in the 



'bulk' of a pure-A fluid, ( [V^^ - V^^] )r, ~ ^ ( [V^^ - V^^] ) 
tion will take the form 

AF(A«) 



A=A(R„) 



AH{X^) 



N 



then the compensa- 



(8) 



Since atomistic and coarse-grained systems usually follow different equations of state 37 
as depicted in Fig. [21 the densities of the two regions will generally differ. It is worth noting, 
though, that by adjusting the number of particles in the system one can easily tune the 
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particle density in the AA region to the reference value p* . In this case the pressure in the 
entire system would adjust to the reference value of the atomistic system. 

A different compensation route has to be taken if, instead of the same pressure, one wants 
to ensure that both subsystems coexist at the same reference density p* . In particular, the 
chemical potential gradient, which is generally established across the transition region, would 
have to be counterbalanced {29 1. This idea leads to the following form of the compensation 



term in Eq. 



A//(A(RJ) . A,(A) = ^ + M^, (9) 

iV p* 



where A/x is the difference in chemical potential across the transition layer and is related 
to the (molar) Gibbs free energy difference by A/i = AG/N = AF/N + Ap/p*. Again, 
Kirkwood TI provides a way to predict Ap by simultaneously evaluating the Helmholtz free 
energy difference AF{X) and the pressure difference Ap{X) in independent simulations of 
pure-A fluids at the reference state (p*, T) and varying A. 

Fig. |2] graphically summarizes the possible routes allowed by these two forms of the free 
energy compensation. The "pressure route", with AH{X) = AF{X)/N, cancels the extra 
interface "pressure" and guarantees that mechanical equilibrium is uniquely established by 
inter-molecular forces (however, in general, the AA and CG subregions will attain different 
densities). On the other hand, in the "density route", the addition of A^r(A) = Ap(A) 
compensates the difference in chemical potential across the transition region leading to an 
equilibrium state where both subsystems coexist at the same density, but different molecular 
pressure. 

To test and validate the proposed compensation schemes we considered the above men- 
tioned tetrahedral system, but we substituted the CG IBl potential with a WCA potential, 
which was deliberately parametrized to give a higher molecular pressure and AF{X) < 
than the atomistic system at the same state point. For sake of simplicity, the same number 
of molecules, total volume of the system and relative volumes of the different subregions 



were used in all simulations (see also |34l|). 

The red curves in Fig. [3] show the pressure and density profiles for the uncorrected H- 
AdResS Hamiltonian, Eq. O Both quantities exhibit jumps in the transition regions, and in 
the AA region neither pressure nor density attain the reference values. Making use of Eqs. 
E] and M we can compensate for the free energy imbalance between the AA and CG regions. 




FIG. 2. Cartoon illustrating the thermodynamics of H-AdResS. The isothermals of AA and CG 
models (black solid and dashed lines) are shown. When no compensation term is added to the 
H-AdResS Hamiltonian in Eq. ^ density and pressure in the two regions are different from 
their reference values (red). Applying compensation terms it is possible to maintain the coupled 
systems either at the same density ("constant density route", blue) or at the same molecular 
pressure ("constant pressure route", green). 

The "constant pressure route" balances on average the drift force (F* ), thus producing a flat 
molecular pressure profile and leading to an average momentum conservation in the whole 
system; the density is nonetheless different in the two regions. In contrast, the "constant 
density route" levels out the density to the reference value p* = N/V by taking the pressure 
in the bulky AA and CG regions to the values they have in the corresponding homogeneous 
simulation. The compensation term of Eq. ([9]) does not take into account density-density 
correlations over the transition layer and, as observed in Fig. 121 this produces small density 
fluctuations (of about ~ 3%) in the transition region. We are currently working on a 
generalization of the present framework to include such correlations. In any case, if required, 
;he small density fluctuations can be removed by an iterative reflnement scheme (see e.g. 



28|). 

To summarize, we have presented a method, H-AdResS, to simulate molecular liquids with 
position-dependent levels of resolution. Whereas in the original AdResS scheme the exact 
enforcement of Newton's 3rd law impedes a general Hamiltonian formulation {20], l^^] , in H- 
AdResS this requirement is relaxed to formulate the problem in terms of a global Hamiltonian 
function. This method allows us to generate equilibrium states in any well deflned statistical 
ensemble, which therefore can be sampled by either Monte Carlo or Molecular Dynamics. 
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FIG. 3. Plots showing the effect of the free energy compensations on the density profile (upper 
panel) and pressure profile (lower panel) in a H-AdResS simulation with CG potential having larger 
molecular pressure than the fully atomistic. The case are the same as to those presented in Fig. 
2. Red, green and blue lines refer to simulations performed using the H-AdResS Hamiltonian in 
Eq. [21 the "constant-pressure" (Eq. [8]) and the "constant-density" (Eq. [9]) compensation routes, 
respectively. All pressures are normalized to the value of the all atom simulation (dash-dot line); 
the dotted line indicates the pressure of the coarse-grained system at the reference atomistic density. 

In H-AdResS, the potential energies of the molecules are weighted according to their 'local 
nature' (atomistic, coarse-grained or hybrid). Based on the analogy with standard Kirkwood 
thermodynamic integration we have proposed two schemes to correct for the drift force 
appearing in the hybrid region. Also these compensation terms are not time- or path- 
dependent, so that no bookkeeping is required to enforce energy conservation; in particular, 
thermodynamical equilibrium is achieved without the help of a local thermostat to remove 
the excess heat produced in the hybrid region. The pressure and density routes for free 
energy compensation offer a simple way to optimize the embedding of the system as well as to 
modulate the thermodynamic balance between AA and CG regions. This new approach thus 
significantly widens the options to couple within a single simulation setup representations 
at rather different resolution, making for a valuable tool for many problems in soft matter 
science. 
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